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Abstract 

We present a molecular dynamics study of force patterns, tensile 
strength and crack formation in a cohesive granular model where the 
particles are subjected to swelling or shrinkage gradients. Non-uniform 
particle size change generates self-equilibrated forces that lead to crack 
initiation as soon as strongest tensile contacts begin to fail. We find that 
the coarse-grained stresses are correctly predicted by an elastic model that 
incorporates particle size change as metric evolution. The tensile strength 
is found to be well below the theoretical strength as a result of inhomoge- 
neous force transmission in granular media. The cracks propagate either 
inward from the edge upon shrinkage and outward from the center upon 
swelling. 

83.80.Fg, 74.80.-g, 45.70.Mg 

key words: granular media / porous media / cohesion / shrinkage / swelling 
/ rupture 



The term "cohesive granular media" covers a vast spectrum of granular ma- 
terials in which rigid grains are bound together by cohesion forces of various 
chemico-physical origins |Maugis, 1999| . Well-known examples are fine pow- 
ders and soils with more or less colloidal or water content |Mitchell, 1993| . The 
solid-like behavior attributed to noncohesive granular media under quasistatic 
shearing |Jaeger et al., 1996| becomes the dominant feature in the presence of 
cohesion, with an increasing effective tensile strength as a function of the con- 
tact tensile strength |Kim and Hwang, 2003| . The stress-stain behavior and 
fracture mechanics of cohesive granular media raise interesting open issues from 
a grain-scale point of view and in interaction with heat or mass transfer. 
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Figure 1: (a) Geometry of the sample; (b) Relative displacements between two 
edge points belonging to two particles and coinciding initially with their contact 
point. 



An appealing issue is how and in which respects these "granular solids" 
differ from molecular solids (in the absence of a granular structure). For 
example, the phenomenon of stress concentration, induced by defects at dif- 
ferent scales, governs the initiation of failure in molecular solids, the effec- 
tive tensile strength remaining generally far below the "theoretical" strength 
|Herrmann and Roux, 1990| . In a granular assembly, stress concentration oc- 
curs already at the particle scale in the form of a highly inhomogeneous distri- 
bution of contact forces |Mueth et al., 1998||Radjai et ai, 1 998 . This suggests 
that, even in the absence of mesoscopic defects, the tensile strength will be 
weak compared to its theoretical value for a granular assembly (to be defined 
below). However, the tensile-strength properties have scarcely been analyzed 
from a microscopic standpoint |Kim and Hwang, 2003| |Radjai et al., 200T| . 

In this Letter, we consider a benchmark test that was designed to probe the 
intrinsic tensile response (reflecting only the granular disorder) of a cohesive 
granular sample by avoiding both wall effects and strain localization as spurious 
sources of randomness. The sample consists of rigid cohesive disks compacted 
numerically into a circular form in a two-dimensional space; See Fig. ^a). At 
the start, the normal force is exactly zero at all contacts. Then, the particle 
diameters are increased (or decreased) at a rate that depends on distance to 
the sample center. Such gradients of particle size change occur, for instance, in 
fine soils, where particle swelling (or shrinkage) happens as a result of humidi- 
fication (or drying) |Mitchell, 1993| . As we shall see in detail below, this bulk 
straining induces a field of radial (or orthoradial) tensile self-stresses increasing 
in magnitude with time, and leading eventually to crack initiation at the center 
(or on the edge). 

For the simulations, we used the molecular dynamics method with 
a velocity- Verlet scheme for the integration of the equations of motion 
|Allen aiid~ Tildeslc y, 1987| . Cohesive interaction between two particles implies 
resistance to relative motion (normal displacement e?„, tangential displacement 
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Figure 2: Tensile (a) and compressive (b) normal forces generated by the 
swelling of a single particle (in black). The line width is proportional to the 
normal force. 



dt and angular displacement 7) of two edge points belonging respectively to the 
two particles and coinciding initially with the contact point; Fig. ^b). The 
corresponding contact actions are the normal force /„, the tangential force ft 
and the contact torque M. We assume a linear elastic behavior characterized by 
three stiffnesses E n , E t and Ey, so that /„ = E n d n , f t = E t d t and M = E y j. 
As usual, damping actions are added in order to account for contact inelasticity 
and ensure numerical stability. 

This elastic behavior holds as long as the contact actions remain below a 
"yield surface" £ = £(/„,/t,M) = 0. We used the following function that fits 
our previous experimental tests with a particular type of glue used to stick 
cylindrical particles together |Uelenne et al., 2002| : 

<-(*)♦(*)' ♦(£)'-*• 

where f%, ff and M v are the yield parameters for normal, tangential and an- 
gular actions, respectively. The elastic domain corresponds to ( < 0. Note 
that /„ can take indefinitely large values (the positive sign corresponding to 
compressive forces) but it has a lower bound /„ = — /* that defines the largest 
tensile force that can be sustained by a contact. As soon as £ > 0, the cohesive 
bond breaks down irreversibly and the contact turns into noncohesive frictional 
behavior |Luding, 1998| . 

In general, the shape of the yield function £ and the values of the parameters 
will influence the failure properties of the material as a whole for a specified 
loading path. In our system, loading by particle swelling or shrinkage induces 
appreciable displacements only along the contact normals. As a result, the 
behavior is not sensitive to the choice of the values of ff and M y , and the 
failure is governed by extensional strain when is reached at a strongly tensile 
contact in the sample. 
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Figure 3: Compressive (a) and tensile (b) forces in a shrinkage simulation. 



We used samples composed of 1133 polydisperse disks with a uniform dis- 
tribution of diameters D within a range [D m i n , D max ] where D max =1.2 D m i n . 
The coefficient of friction is \i = 0.1. Each sample is created by removing all 
particles belonging to a noncohcsive assembly in static equilibrium except those 
contained inside a circle of radius Rq . The cohesive bonds are then switched on 
and the sample is allowed to relax to equilibrium. The samples prepared by this 
procedure correspond to a dense but disordered packing of solid fraction ~ 0.89 
and coordination number 3.8. 

At the beginning, the system is in its reference state with the contact actions 
being identically zero everywhere. Obviously, multiplying all particle diameters 
by the same factor does not disturb this state since no relative motions are 
generated at the contact points, whereas differential particle-size change gives 
immediately rise to permanent compressive and tensile force gradients. For 
instance, the swelling of a single particle produces compressive radial forces by 
pushing the neighboring particles outward, as well as tensile orthoradial forces 
by increasing the total length of the "rings" of contiguous particles surrounding 
the swelling particle; Fig. [21 A slight shrinkage of the same particle produces 
exactly the same force patterns with the signs inverted everywhere (compressive 
contacts turning to tensile, and vice versa). 

Since we are interested here only in the effect of bulk straining, we require 
that the swelling rate Si = Di/ Di of each particle i is independent of its diameter 
Di. We use the simplest swelling kinetics defined by a constant gradient from 
the center to the edge, Si = fj, where r, is the distance of the particle to 
the system center, and a is a constant rate. Positive and negative values of a 
correspond to particle swelling and shrinkage, respectively. 

Figure shows snapshots of normal compressive and tensile forces in a 
shrinkage simulation. Although at the very local scale the forces are inho- 
mogeneously distributed, we observe radial and orthoradial compressive forces 
decreasing in magnitude from the center to the edge, as well as orthoradial 
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Figure 4: Crack patterns in swelling (a) and shrinkage (b) simulations. 




tensile forces decreasing in magnitude from the edge to the center. The cracks 
appear on the edge as soon as the first tensile contact fails, and they propagate 
towards the center as shown in Fig. 0Jb). In swelling simulations, the respective 
roles of compressive and tensile roles are simply interchanged with respect to 
the shrinkage case. As a result, the cracks are initiated at the center, and they 
propagate towards the edge; Fig. 01 a). 

The coarse-grained stresses can be evaluated from the contact forces 
by means of the "micromechanical" expression of the stress tensor cr 
|Christoffersen et al., 1981] : 

= h E st 3> (2) 

c£V 

where i and j design the coordinates, V is the control volume over which the 
stress tensor is evaluated (the contacts c taken from this volume), is the i 
component of the force f c at contact c, and £j is the j component of the vector 
t c joining the centers of the two contact neighbors. 
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Figure 6: Evolution of the sample radius R in swelling (a) and shrinkage (b) 
simulations. 



The stress tensor is a well-defined average if the control volume V contains 
a sufficiently large number of contacts. This requirement is satisfied by taking 
concentric circular volume elements as suggested by the rotational invariance of 
our system. We use polar coordinates and the radial and angular positions will 
be denoted by r and 8, respectively; Fig. QJa). As a result of isotropic straining, 
the cross components o~ r g are zero. 

The radial stress a rr and orthoradial stress o~gg are shown in Fig. [5fa) as a 
function of distance r to the center at different stages of a swelling simulation. 
For each data set at a given instant of evolution, we have normalized the distance 
r by R and the stress components by the largest tensile stress — cr max (occurring 
at the center). We see that the normalized stresses collapse on a straight line 
as a function of r and they nicely agree with a one-parameter analytical fit that 
will be detailed below. Note that a rr is negative (tensile stress) throughout 
the sample, whereas age changes sign at r ~ R/2. In the case of shrinkage 
simulations, similar results are obtained with opposite signs, as shown in Fig. 

We now turn to analytical evaluation of the stresses. At a coarse-grained 
scale, the granular assembly will be represented by a linear elastic medium with 
an effective stiffness E and an effective Poisson's ratio v. This is a plausible 
assumption although, as we shall see below, the behavior of the stresses as a 
function of r is independent of the nature of the interactions. On the other 
hand, coarse-grained swelling at a point A of polar coordinates (r, 9) in space is 
equivalent to an imposed isotropic space dilation s(r, 6) — (si)i e v{r,e) > where the 
average is taken over all particles contained in a representative volume V{r, 6) 
centered on A. Since the rates are independent of particle diameters, we get 



Then, the strain-rate tensor e at a point is the sum of two terms: 

i = e e + si, (4) 
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where i e is the elastic strain rate, J represents the unit tensor, and si is the 
metric change rate. 

We assume that the displacement field u(r, 8) is radial (according to the 
symmetry of straining expressed by Eq. and that of the sample) so that 
e r g = 0. Since the radius R of the sample changes with time, Hooke's laws will 
be written in the form of rate equations: 



where extensional strains and compressive stresses are counted positive. We 
need also the balance equation which takes the following form in polar coordi- 
nates: 



The set of equations and is easily integrated over time and space using 
the boundary conditions u(r = 0) = (imposed byOJ) and a rr (r — R) = (by 
continuity of the normal stress at the boundary). The solution is 



We see that both stress components are linear in r. The simulation data of Fig. 
El were fitted by adjusting only the effective elastic modulus E. The evolution 
of the system is, however, nonlinear as a function of time. The evolution of 
R is shown in Fig. for swelling (a > 0) and shrinkage (a < 0) simulations 
together with the analytical fit from Eq. [7| which involves no fitting parameter. 
The agreement is excellent although the nonlinear nature of the evolution can 
not be seen for \a\t <C 1. The largest tensile stress a max occurs on the edge for 
shrinkage and at the center for swelling. From Eq. [7| we get a max = |-E \a\t. 
Again, this linear form nicely fits the evolution of a max (by virtue of the fits 
shown in Fig. |3J) up to failure for a max = a v . The latter represents the effective 
tensile strength of the material. 

It is worth noting that, Posisson's ratio v does not appear in Eqs. 0and 
the only role of the stiffness E is to set the stress scale. This means that the 
behavior of the stress components and and sample size as a function of r is 
independent of the local force law. In particular, in the limit of infinitely rigid 
particles, the same results remain true up to a stress scale which may be fixed 
through a confining pressure. More generally, both the local interactions and 
the mass or heat transfer influence the stress scale. 

By analogy with molecular solids, we introduce a "theoretical" 
tensile strength u y th based on the interactions between two particles 
|Herrmann and Roux, 1990| . According to Eq. [21 the orthoradial stress is 
ogg = n c (fgig) — n c (£) (fe), where n c is the density of contacts and (. . .) 
designs averaging over the control volume. The largest value of agg in tension 
corresponds to the limit where all forces are polarized in the same direction and 



E% r = S rr — S = —^(& rr — l/&gg), 
igg = Egg — S = — ^ (<Tgg — V & rr ), 



(5) 





7 



they have all reached the largest tensile force fv. This defines the "theoretical" 
tensile strength 

<4 = MQfZ (8) 

In our simulations, the measured tensile strength a y is by a factor ~ 4.3 
below <j\ h . In molecular solids, a similar discrepancy between of h , defined 
from atomic interactions in a regular atomic arrangement, and o~ y stems 
from "built-in" disorder at different scales leading to stress concentration. 
In a granular solid, the disorder is "intrinsic" to the structure. As a re- 
sult, the contact forces both in cohesive and noncohesive granular media 
have a wide distribution with a decreasing exponential shape for strong forces 
|Radjai et al., 1998| |Radjai et al., 2001] . It seems that a y reflects the strongest 
contact force at failure, whereas a v th represents the mean tensile force by con- 
struction. Indeed, the strongest tensile force in simulations is by a factor ~ 5 
below the mean tensile force, and this is close to a y h /a y ~ 4.3. As far as we 
know, this is the first example showing how the local force inhomogeneities in a 
granular material control a macroscopic property, namely the tensile strength. 

In summary, our numerical data and their comparison with an analytical 
evaluation of stresses in the elastic domain and at failure suggest that a macro- 
scopically elastic behavior is relevant up to crack initiation, as in molecular 
solids. However, the tensile strength is dependent on the inhomogeneous trans- 
mission of forces. The simple test described in this Letter not only provides 
reproducible results, but it has also the advantage of combining features of 
discrete modeling with theoretical predictability at the macroscopic scale. 

This approach may now be used to investigate and to predict the tensile 
thresholds and crack propagation in cohesive granular materials as a function of 
the intial density and anisotropy of the material or the possible couplings of the 
local cohesion with mass and heat transfer in the pores as in fine soils and gran- 
ular rocks |Tarbuck and Lutgens, 2002] , The theoretical approach can be ex- 
tended to other structured media involving mesoscopic length scales, such as gels 
|Mrani et al., 1995] , cellular media |Schwarz and^afran^2j)02j, layered struc- 
tures such as wood |Kubler, 1987| and pastes |Ponsart^t^bT^2003| . Swelling or 
shrinkage may occur as a result of cellular growth (in biological systems) or the 
evolution of local variables such as water content and temperature. 

It is a pleasure to thank J.-C. Benet and J. N. Roux for helpful suggestions. 
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